Indirect Quantum Tomography of Quadratic Hamiltonians 
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A number of many-body problems can be formulated using Hamiltonians that are quadratic in 
tiie creation and annihilation operators. Here, we show how such quadratic Hamiltonians can be 
efficiently estimated indirectly, employing very few resources. We find that almost all properties of 
the Hamiltonian are determined by its surface, and that these properties can be measured even if 
the system can only be initialised to a mixed state. Therefore our method can be applied to various 
physical models, with important examples including coupled nano-mechanical oscillators, hopping 
fermions in optical lattices, and transverse Ising chains. 



I. INTRODUCTION 



There has been considerable interest in the problem of Hamiltonian identification through indirect probing , thereby 
developing various quantum mechanical versions of classical system tomography or classical 'inverse scattering' prob- 
lems For certain types of interactions, it was found 0-131 that only few resources are required to obtain an accurate 
model of the system. Indirect Hamiltonian estimation is therefore an interesting problem for both pragmatic purposes 
and fundamental insights. We are interested in the following questions: How can we obtain precise information about 

■ a Hamiltonian under restricted access? What can we learn about the 'inside' of a large system by only looking at a 
subsystem of it? Under which conditions is such indirect probing possible? 

Recent studies have focused on this problem for cases of chains and networks of spin-1/2 particles. The common 
question addressed can be formulated as follows: Can we estimate all parameters, such as coupling strengths and 
^2 ' local fields, by accessing only one or a few spins? It should be emphasised that even direct Hamiltonian estimation, 
or more generally, process tomography, is hard, because the required number of measurements and the complexity of 
the post-processing both scale exponentially with the system size. However, in realistic situations, we usually have 
I ] a priori knowledge based on the underlying physics. It has been shown that such knowledge can be used to develop 

■ compressed sensing protocols [1, 0| , which greatly reduce the complexity of process tomography. Various works on 
' indirect Hamiltonian estimation have relied on similar assumptions; namely, that the dynamics is restricted to a 

T-H subspace of polynomial dimension [2-4]. In [2J, the efficiency of the estimation in terms of the required time and the 
number of measurements is discussed. An interesting example that does not rely on a subspace was analysed by Di 
, Franco et al. We will see here that this is a special case of the generic estimation of quadratic Hamiltonians, 
• which can be estimated efficiently due to a simple description of their dynamics in the Heisenberg picture. Di Franco 
also found that the estimation is quite robust against noise In Q the ID methods were generalised to arbitrary 
graphs, and the possible elimination of degeneracies was discussed. Also, Wiesniak and Markiewicz |4j] went beyond 
I ' the simplest subspace in order to study quasi-lD systems. Table |T] summarises the results obtained so far in terms of 
^ the settings and assumptions considered. However, the analysis of physically important cases, such as the transverse 
Ising model and the XY model with a magnetic field, have remained open. Solutions to both cases will be presented 
... in this work. 

' Our main goal in this paper is to develop a method to perform indirect quantum tomography for many-body systems 
of identical particles. Even though the method is analogous to the spin case, the Hamiltonians considered here have 
a higher number of parameters, and it is surprising that they can still be estimated in a similar fashion. The class of 
Hamiltonians we study here is those of quadratic form in bosonic or fermionic operators. There has been tremendous 
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Figure 1: Indirect classical system tomography of a quadratic Hamiltonian. In this example, the spring constants ki and masses 
nii of a chain of coupled harmonic oscillators can be determined by monitoring the dynamics of a single particle at the end (in 
red). See [l| for details. 
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Interaction type 


Needs preparation 


Geometry 


Obtain 


Reference 


XX + YY + AZZ 


specific state 


ID 


couplings 


(2] 


(l + 7)XX + (l-7)Kr 
7/1,-1 


no 


ID 


couplings 


[5] 


XX + YY + AZZ + Z 


specific state 


arbitrary 


couplings and fields 


(31 


XX +YY+Z 


specific state 


quasi-lD 


couplings, partial topology 


m 


a'' a + aa + h.c. 
(fermions or bosons) 
{l + j)XX + {l-j)YY + Z 


arbitrary state 


arbitrary 


couplings, fields, anisotropics 


this paper 



Table I: Overview of indirect Hamiltonian estimation schemes. The interaction types represent the Pauli matrices involved in 
the spin-spin coupling, e.g., XX + YY + AZZ + Z stands for a Hamiltonian of the form ^ Anm {XX + YY + AZZ)^ ^ + 

progress in experiments of quantum random walks [l^l, optical lattices [ll|, coupled cavities nano-mechanical 
oscillators etc., which can be modelled by such quadratic Hamiltonians. Thus, the indirect estimation scheme 
we present here will be of use in reducing the necessary resources for modelling such systems. For the case of bosons, 
it is the most direct translation of the work by Gladwell [H to the quantum case. Gladwell studied how the spring 
constants and masses of coupled classical harmonic oscillator chains can be estimated by looking at the movement 
of only one particle (see also Fig. [1]). Furthermore, our estimation protocol gives a natural generalisation of the 
problems on spin chains, since quadratic Hamiltonians of fermions also describe a certain class of spin systems, such 
as the transverse Ising model. 

The main method of indirect estimation is summarised as follows. First, the system is initialised to an arbitrary 
but fixed state. This can be, for example, even a thermal state, and can occur on slow time-scales via relaxation. 
Then, some simple single-particle properties are initialised and measured at a later time. Finally, the accumulated 
data is Fourier transformed, and the parameters are extracted through a set of linear equations. This simple method 
is outlined in Fig. [2] for the ID case. The procedure is analogous to an 'inverse scattering' problem because the 
perturbation introduced in one edge of the sample (e.g., rotation of the first qubit) propagates through the sample, 
'scattering' with the inner structure of the Hamiltonian, and then this information encodes the structure of the system. 
While it is obvious that this procedure provides some information on the system, the surprising result here is that 
all information can be uniquely identified. Since we obtain information on anisotropics as well, not only the coupling 
strengths but also the type of interaction can be determined by our method. 

The paper is structured as follows. First, we introduce the necessary notation for quadratic Hamiltonians and some 
techniques for their diagonalisations. Since these are well established methods, we will only discuss them as a 'recipe' 
for the estimation procedure. We then discuss the simplest case of estimation, namely when the system is a chain of 
hopping particles, and later generalise it to arbitrary graphs. Finally we discuss how the results apply to ID chains 
of spins and conclude. 

II. NOTATION AND DIAGONALISATION 

The most general quadratic Hamiltonian of N indistinguishable particles is written as 

JV 1 ^ 

n,m— 1 n,m— 1 

where the a]j and a„i are creation and annihilation operators and A, B are matrices describing the parameters we 
would like to estimate. For H to be Hermitian we must have A = and = —eB, where we introduced the 
parameter e = 1 for fermions and e = — 1 for bosons. We will mostly follow the notation of though we shall write 
all vectors in Dirac notation. At first, we put all Hamiltonian parameters into the Hermitian 2N x 2N matrix 



(2) 
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Fourier transform to obtain Ek and {l\Ek) 



Solve linear equations to obtain H 



Figure 2: Schematic overview of our estimation sciieme for a ciiain of lengtli A*'. Tiie main requirements are (i) relatively fast 
single-qubit rotations and measurements on the first site, and (ii) a decoherence time of at least 7V^. 



and introduce the column vector operator 



so that Eq. ^ can be expressed up to a constant as 

H = {l/2)a^Ma. 

Throughout this paper, we make two technical assumptions: first, all the phases of M^- are assumed to be known. 
Although some phases might be easy to determine, this requires complicated studies of gauge invariance that do not 
seem to be worthwhile, as in many practical cases all elements are real and positive. Second, for the bosonic case 
we assume that the matrix M is positive definite. Again, in principle this can be generalised, but this way we avoid 
technical difficulties of symplectic transformations (l4| . 

As in the efficiency of our method depends on how many entries of M are a priori known to be zero; that is, 
on knowledge of the coupling graph. If such knowledge is not available, we have to perform measurements on all but 
one of the qubits. If the graph is known to be highly sparse (for instance a chain) we only need to access a single 
qubit. But before going into the details of Hamiltonian identification, let us briefly review the diagonalisation of 
the Hamiltonian Eq. ([1]) and thus the dynamics, introducing some notations. For more detailed descriptions on the 
diagonalisation procedure, see, e.g., (T^ . 

For quadratic Hamiltonians of the form in Eq. ([T]) , there exist quasi-particle creation and annihilation operators foj. 
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and bi, with which the Hamiltonian can be represented by the simple form of non- interacting modes, 

H = ^ Ekbl,bk + const. 



(3) 



For this reason, quadratic Hamiltonians are also referred to as 'quasi-free' interactions. We need to know the trans- 
formation T that maps the operators a and for particles to b and b^ for quasi-particles, i.e., /3 = Ta, where f3 is 
defined by 



In order to ensure the canonical commutation relations for the operators bk and bj., T must satisfy = rjT^rj, where 

I el i 



The Hamiltonian is now written as 



H 



-f3''r]{TT]MT-^)(3 + const. 



It can then be shown that rjM is diagonalised by T as 



Tr]MT~^ = 



E 

-E 



where E — diag {Ei , E^}, to have the desired form of Eq. ([3]). Note that the energy eigenvalues appear in pairs of 
positive Ek and negative E^+n = ~Ek values (fc = 1, ...,N). 
The matrix T consists of the right eigenvectors \Ek) of rjM as 



/ {Ei\ 



T = 



\{E: 



2N\ 



and the inverse of T is given by 



\Ei) 



\E2n)). 



For bosons, the matrix rjM is not Hermitian and the distinction between right and left eigenvectors is necessary. This 
gives rise to a few further peculiarities, such as the modified normalisation and completeness relationship (see below). 
Fortunately, in this work we are solving an inverse problem, and do not have to discuss how to find these vectors and 
how numerically stable the corresponding algorithms are. 

It is worth pointing out that the \Ek) are not representing physical states but just introduced here as a part of 
solving the Heisenberg equation of motion for the creation and annihilation operators. The completeness relationship 
is given by 



N 



J2\Ek){Ek\v + ^\Ek+N){Ek+N\v ^ ^ 



2Nx2N, 



(4) 



fe=i 



and the vectors \Ek) are chosen to fulfil the normalisation relationship 



{Ek\v\Ek') = Vkk' 



(5) 
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For convenience, let us also introduce vectors |n) as the canonical basis vectors: 



|n) = 1 -f- n-th row. 


Vo/ 

Due to the structure of the matrix M, the upper and lower eigenvectors of rjAI are related as 

where © is the addition modulo 2N. The dynamics of the original operators a can be found from /3„(t) = e~*^'"*/3„(0) 
as 



(6) 



where we have introduced a sign function s(m, fc) through 



s{m,k) = 1 (m = l,...iV;fc = 1,...7V) 



s(to, k) 
3(111, k) 



{m = l,...N;k^ N + l,...2N) 
{m = N + l,...2N\k^l,...N) 



s{m,k) = 1 {m = N + l,...2N]k^ N + l,...2N). 



III. ESTIMATION OF CHAINS 



A. Experimental requirements 

Let us first consider a ID chain of interacting particles, meaning that A and B are tridiagonal. Also assume that 
we can initialise the chain in a fixed state po- This state could, for instance, be a thermal state, but we do not require 
to have the exact form of p^: we just have to be able to repeatedly initialise the chain to the same state po- As 
before we perform initialisations followed by measurements at the first site. The quantity we need to measure 

at the first site is oi ± a\ for different times up to iV^ In order to eradicate the dependence on the initial state 
we measure two sequences of (ai(t)) after preparing the first site to give two different initial values, i.e., (ai(0)) — ci 

and (ai(0)) = C2. Using (a„) — (^a[^, and thus (a|(0)) = c*. and subtracting the measurement results, we obtain a 
quantity that only depends on Ac = ci — C2 . It is given through Eq. ([6]) by 



27V 



.(0)) 



m,fc— 1 
2N 





2N 




E - 


- 1 


m,fc— 1 



27V 



Ac 5] .(1, fc)e-^'=*|rrfe'P + Ac* J2 <N + 1' k)e-^'^''%-k (T^"')^^! ■ 



(7) 



fe=i 



fe=i 



This initialisation can be performed by a von Neumann measurement, or, as long as the reduced density matrix at 
site one is not maximally mixed, by applying different single qubit rotations (in some experiments von Neumann 
measurements are hard). As we see in Eq. the dependence on the initial state is completely removed. This is 
thanks to the absence of the interactions between particles: they (almost) do not see each other, so the information 
on the 'injected' particle can be extracted by subtracting the influence from others. 

For the spin chain case the eigenfrequencies are non-degenerate and T^j} = 7^ (Vfc) 0, @|. For the present 

case of quadratic Hamiltonians we were unable to prove this, but could confirm it numerically. Hence, a Fourier 
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1) |2) |3) |4) \N-1) \N) 




\N + 1) \N + 2) \N + 3) \N + 4:) \2N - 1) \2N) 

Figure 3: The matrix M as a coupling graph on the nodes \n). 

analysis provides us with the frequencies Ek and the amphtudes Ac |TYfe^p + Ac*eT^f} (T^^)^ . Summing these 
amphtudes gives the value of Ac : 

N 2N 

$:(Ac|T-T + Ac*eT-(r-);^J+. (^c I^^P + Ac*eT- (T^);^^^) = 

fe=l fc=A'+l 

Ac(l|l) + Ac*e(l|iV + 1) = Ac, 

where we used the completeness relationship Eq. Equation ([7]) still contains mixtures of the coefficients IT^^I and 

"^ik N+i ■ '^^^ separate them by measuring another pair of initialisations c[ : as long as Ac' ^ rAc (r S M), 

we can solve the linear equation for iT^Jj^^j. Without loss of generality, we choose |(l|i?fe)| = l^'ifc^l — T^j^ — (l|i?fe) (Vfc) 
by arranging the global phase of each eigenstate \Ek)- In conclusion, a few random rotations or initialisations of the 
first qubit, followed by measurements, provide us with E}. and {l\Ek). 

Let us now describe how to obtain the parameters of M from this observed data. We have to distinguish between 
the generic case where the off-diagonal couplings An^n+i and B„^n+i are distinct, and the special case where they are 
equal. 

B. Different off-diagonal couplings 

As we have seen above, what we diagonalised is the 2N x 2N matrix rjM, so it is helpful to regard M as a 
representation of a graph consisting of 2N nodes (see Fig. [3]). Its off-diagonal entries correspond to the coupling 
strengths between nodes, whereas the diagonal elements represent the intensity of the 'field' at each node. We can 
then start with a recursive algorithm similar to [l|-[l| by applying M to the localised states at sites 1 and + 1; 

M\l) = An\l} + A2i\2) - eB;,\N + 1) - eB;,\N + 2) 
M|iV + l) = Bn\l)+B2i\2)-eAl,\N + l)-eA;,\N + 2). 

Using these equations for {l\riM\Ek) and {N + l\riM\Ek) we arrive at 

Ek{l\Ek) = {l\vM\Ek) = {l\M\Ek) - Aii{l\Ek) + A;^{2\Ek) - eBn{N + l\Ek) - eB^iiN + 2\Ek) 
eEk{N + l\Ek) = e{N + l\r^M\Ek) ^{N + l\M\Ek) = B*^{l\Ek) + B;^{2\Ek) ~ eAn{N + l\Ek) - eA^iiN + 2\Ek) . 

Note that Bn = for fermions. Among these parameters, we can immediately obtain ^nand Bn from the quantities 
that are estimated by measuring {ai{t)) as 

N 

An = {l\M\l}=Y,[Ek\{l\Ek)\'' + eEk+N\{l\Ek+N)f], 

k=l 

N 

Bn = {l\vM\N + 1) = - ^ [Ek{l\Ek){Ek\N + 1) + eEk+N{l\Ek+N){Ek+N\N + 1)] . 
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Therefore we can collect all the known terms of the above equations on the left-hand-side as 

known = Al^{2\Ek) - eB2i{N + 2\Ek) (8) 
known = Bl^{2\Ek) - eA2i{N + 2\Ek). (9) 

Taking linear combinations, we obtain 

known = 
known = 

It may appear as if the right-hand-side contains too many unknowns to solve these equations. However, similar to 
the original work by Gladwell [l| and the spin chain case, 0|, we can use the normalisation of the eigenstates. In 
this case, it is given by Eq. ([5]). Summing up the mod squares of the above equations, the dependencies on {2\Ek) 

and {N -\- 2\Ek) vanishes and we can infer the absolute value of each coefficient, i.e., ^|^2i|^ — |-B2i|^^ and 

\ST2\ (l^2i|^ ~ |-S2i|^y If |^2i| 7^ \B2i\, we obtain |i?2i/^2i| = 3 by dividing the above two equations and then 

through 1^21 1(1 — g^), both IA21I and |i?2i| are obtained. Because the phases of A and B are known, then we learn 
A21 and 321- We can then express a similar set of equations for the next site {2\rjM\Ek) and {N + 2\rjM\Ek) ■ By 
induction, we then see that all matrix elements of A and B can be obtained, as desired. It is worth pointing out that 
the scheme works even though the graph in Fig. [3] is not infecting Q. 

We will now look at the cases with equal off-diagonal couplings in more detail, because such physical systems are 
often encountered, e.g., transverse Ising for fermions, coupled harmonic oscillators for bosons. 



C. Equal ofF-diagonal couplings 

When |A„+i.„| = |_B„+i^„|, the above method fails. This is the case for interacting harmonic oscillators without the 
rotating wave approximation [l5| and for quantum Ising models, and therefore of interest in a number of practical 
situations. We focus on the case where A and B are real, e.g., An^n+i = ^n.n+i = —^Bn+i.n- The diagonal elements 
Ann and _B„„ are always different if there is a transverse field (fermions) or if the masses are finite (bosons), so it is 
reasonable to assume A„„ 7^ i?„„ (the Ising model without transverse field cannot be estimated using our method 
because excitations do not propagate). 

It is convenient to introduce 

|n±) = i=(|n)±|7i + iV)),n = l,...,^. 

A simple calculation then shows that 

Mr;|n±) = (1 ± e)A„,„_i|n - 1^) + (1 T l)A.+i.nbi + 1^) + [Ann ± S„„) jn^), n = l,...,N, 

where we set Aqi = (^v-i-i) = 0- In some sense, this is similar to a ID chain case. As the elements are 
already known, we learn An and Bn from 

(lT|r7M|l±) = ^ii±i-^Bii, 

which can also be written in terms of the known variables by inserting the completeness relation Eq. (|4]). Using 
Ek{l^\Ek) = {l^\rjM\Ek) and normalisation we obtain A21 and |2+). For bosons, we then obtain A22 — B22 = 
{2'^\rjM\2~) through the completeness relation, followed by {2~\Ek) through 

Ek{n+\Ek) - {n+\T^M\Ek) = (A„„ - Bnn){n-\Ek). (10) 

and finally 



A22+B22 = (2-|?7M|2+), 
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by completeness again. On the other hand, for fermions, 

Ek{2+\Eu) = {2+\riM\Ek) = 2A2i{l-\Ek) + A22{2~\Eu) , 

and through normahsation we obtain A22 and {2^\Ek). Knowing all parameters at site 2, we can then proceed through 
induction. For the most general case, we can also allow for chains which sometimes have equal off-diagonal couplings 
and sometimes different ones, by alternating between the strategies described here and in the last subsection. 



D. Estimation of general graphs 

We now briefly describe how the linear case is generalised to arbitrary graphs. This is almost identical to [sj, so we 
will not repeat the details. Similar to the spin case, in the general graph setting, measurements on a single spin do 
not suffice: we need to consider transport in the network. Depending on the network topology, we choose a set C of 
'infecting' [s^ nodes, which are the ones we will perform initialisations and measurements on. For clarity, let us recall 
the definition of graph 'infection'. Suppose that a subset C of nodes of the graph is 'infected' with some property, e.g., 
the flu. This property then spreads, infecting other nodes, by the following rule: an infected node infects a 'healthy' 
(uninfected) neighbour if and only if it is its unique healthy neighbour. If eventually all nodes are infected, the initial 
set C is called 'infecting'. 

Similar to the measurements described in Subsection IIII Al initialising the site m e C and measuring the £th node 
after some time, we can obtain 

{a,{t)) = ^.(m,fc)e-^'=*r,-.i(T,;i)t(a„(0)) 

■m.k 

k k 

= s{m, fc)e-^'=*r,-.i (T,;i) t (a„(0)) + ^ ^(m, fc)e-^'*r,-.i (T,-!) t (at„(0)) . (11) 

k k 

Again, the dependence on the initial state po may be removed by subtracting data for different initial conditions on 
the sites m and i. Starting from some element in C, say m = £ = 1, we can get T^f} as described in Subsection IIII Al 
Then we initialise in m and measure at a different site £ £ C, obtaining including its phase from Eq. pip . Hence, 
all Tki with k^i € C can be learnt from simple experiments on the set C. The ID estimation and infection are then 
used to infer the remaining parameters, as described in more detail in 



IV. APPLICATION TO ID SPIN CHAINS 



Naturally, the above scheme can be applied directly to many cases of Hamiltonian identification for systems of 
spin-1/2 particles. A typical example is the XY chain of spin-1/2 particles, 

N N 

H = Y, c„,„+i[(i + i)s:s:+, + (1 - 7)5^5^+1] + E ^-^n, 

n—l n—1 

as it can be transformed into quasi-free fermions by means of the Jordan- Wigner transformation. As it has been 
noted already in the estimation for this model can be done without initialisation of the entire chain. In the 
Jordan- Wigner picture this becomes very clear. That is, after locally measuring an eigenstate of Xi = ai + a[, thus 
making (Zi) = 0, the initial expectation values of the a„ and aj^ n > 1 are all zero. This is because the Jordan- Wigner 
transformation of a„ {n > 1), i.e., a„ = (^nY\m<n^^^ always contains Zi in the product. One might say that the 
local initialisation in the spin picture corresponds to a global initialisation in the fermionic picture. Combined with 
the weak dependence of local observables on the initial condition that comes from the quasi-free interaction, the state 
dependence of the measurements at the first site is completely removed. Hence our scheme is a proper generalisation 
of (5^1 to include magnetic field and the transverse Ising case, which is important in various physical systems, e.g., 
superconducting (fiux) qubits [16.1 , NMR, etc. Such models have also attracted attention in the context of indirect 
quantum control recently (itI. Il8|. where our estimation scheme is crucial. 
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V. CONCLUSIONS 

We found a simple and efficient method to identify the Hamiltonian of a system of coupled bosons or fermions. While 
the methods are completely analogous to the spin case it is surprising that the higher number of parameters in 

the Hamiltonian that arises from the non-conservation of excitations can still be estimated using the same resources. 
Similarly to we can deal with very weak system initialisation, such as thermal states. Therefore our methods can 
drastically reduce the required resources for system identification. Since we allow for site-dependent anisotropics, not 
only the coupling strengths but also the type of interaction is determined along the way. From the theory side, our 
work once more confirmed a type of 'area law' for estimation: looking only at the surface of short-range interacting 
systems can completely determine their Hamiltonian. It would be intriguing to see if this has direct connections with 
area laws of entanglement [l^. While the efficiency of our method relied on the quadratic form of the Hamiltonian, 
we conjecture that even for models with true interaction terms, e.g., quartic terms in the Hamiltonian, all system 
parameters remain discoverable on the surface in principle. 
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